Tarea 1: Foundations

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en análisis exploratorio de datos y conceptos introductorios de probabilidades. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de ejercicios prácticos con el fín de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Elaboración de código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Pregunta 1: Visualización de Datos

Para esta pregunta usted deberá trabajar en base al conjunto de datos hearth_database.csv, el cual esta compuesto por las siguientes variables:

En base al dataset propuesto realice un análisis exploratorio de los datos (EDA). Para su análisis enfoquen el desarrollo en las siguientes tareas:

Respuesta

# Importamos dataframe
my.frame1 <- read.table(file="hearth_database.csv", header=T, sep=",")
my.frame <- na.omit(my.frame1)

# Vemos nombres de columnas
columns <- colnames(my.frame)
print(columns)
##  [1] "target"   "sex"      "fbs"      "exang"    "cp"       "restecg" 
##  [7] "slope"    "ca"       "thal"     "age"      "trestbps" "chol"    
## [13] "thalach"  "oldpeak"
# Para ver tipos de los datos
type.columns <- sapply(my.frame, class)
print(type.columns)
##      target         sex         fbs       exang          cp     restecg 
## "character" "character" "character" "character" "character" "character" 
##       slope          ca        thal         age    trestbps        chol 
##   "integer"   "integer"   "integer"   "integer"   "integer"   "integer" 
##     thalach     oldpeak 
##   "integer"   "numeric"
# Para clasificarlos según nuestras necesidades
classify <- function(column) {
  if (is.numeric(column)) {
    if (length(unique(column)) == 2) {
      return("Binaria")
    } else{
      return("Numérica")
    }
  } 
  else if (is.character(column)) {
    return("Categórica")
  } 
  else if (is.logical(column)) {
    return("Binaria")
  } 
  else {
    return("Otro")
  }
}

# Aplicamos classify al df
type2.columns <- sapply(my.frame, classify)
print(type2.columns)
##       target          sex          fbs        exang           cp      restecg 
## "Categórica" "Categórica" "Categórica" "Categórica" "Categórica" "Categórica" 
##        slope           ca         thal          age     trestbps         chol 
##   "Numérica"   "Numérica"   "Numérica"   "Numérica"   "Numérica"   "Numérica" 
##      thalach      oldpeak 
##   "Numérica"   "Numérica"
# Tomamos las columnas que son del tipo numérico
columns.use <- names(type2.columns)[type2.columns == "Numérica"]
print(columns.use)
## [1] "slope"    "ca"       "thal"     "age"      "trestbps" "chol"     "thalach" 
## [8] "oldpeak"

En la casilla anterior se clasifican las columnas por tipos “binario”, “categórico”, “numérico”, “otro”. Se decidió que las columnas que tuvieran valores numéricos “limitados” (por ejemplo que toman valores 0, 1 o 2) igualmente se considerarían numéricos (por ejempo “ca”).

# Primera parte
means <- sapply(columns.use, function(col) mean(my.frame[[col]]))
print(means)
##       slope          ca        thal         age    trestbps        chol 
##   1.3993399   0.7293729   2.3135314  54.3663366 131.6237624 246.2640264 
##     thalach     oldpeak 
## 149.6468647   1.0396040
medians <- sapply(columns.use, function(col) median(my.frame[[col]]))
print(medians)
##    slope       ca     thal      age trestbps     chol  thalach  oldpeak 
##      1.0      0.0      2.0     55.0    130.0    240.0    153.0      0.8
quantiles <- sapply(columns.use, function(col) quantile(my.frame[[col]]))
print(quantiles)
##      slope ca thal  age trestbps  chol thalach oldpeak
## 0%       0  0    0 29.0       94 126.0    71.0     0.0
## 25%      1  0    2 47.5      120 211.0   133.5     0.0
## 50%      1  0    2 55.0      130 240.0   153.0     0.8
## 75%      2  1    3 61.0      140 274.5   166.0     1.6
## 100%     2  4    3 77.0      200 564.0   202.0     6.2
maxs <- sapply(columns.use, function(col) max(my.frame[[col]]))
print(maxs)
##    slope       ca     thal      age trestbps     chol  thalach  oldpeak 
##      2.0      4.0      3.0     77.0    200.0    564.0    202.0      6.2

Para la primera parte se aplicó la función sapply para calcular los valores pedidos. Se usó el resultado obtenido del análisis anterior para determinar a qué columnas aplicar las métricas.

# Segunda parte
corr.matrix <- cor(my.frame[columns.use])
heatmap(x = as.matrix(corr.matrix), Colv = NA, Rowv = NA)

Aquí se tiene la matriz de correlación de pearson entre las variables. Se observa la diagonal con valores iguales a 1(comparación de la variable consigo misma) y luego, los valores entre el resto que varía entre -1 y 1.

# Tercera parte
boxplot(my.frame[columns.use[4]],
        main = "Boxplot for Ages",
        ylab = "Ages",
        names = columns.use[4]
)

boxplot(my.frame[columns.use[5:7]],
        main = "Boxplots for Numeric Attributes",
        xlab = "Attributes",
        ylab = "Values",
        names = columns.use[5:7]
)

boxplot(my.frame[columns.use[1:3]],
        main = "Boxplots for Numeric Attributes",
        xlab = "Attributes",
        ylab = "Values",
        names = columns.use[1:3]
)

boxplot(my.frame[columns.use[8]],
        main = "Boxplot for oldpeak",
        ylab = "Values",
        names = columns.use[8]
)

En esta parte se fueron separando los boxplots según los rangos de valores que abordaban los atributos para no afectar la correcta visualización de esto a causa de la escala. Para el caso de las variables que son numéricas pero con un muy limitado rango de números naturales (slope, ca y thal) hay outliers, sin embargo, deben ser valores que no se repiten muchas veces y que sí son válidos.

# Cuarta parte
# Hay que dividir el dataset según el target y hacer los histogramas
T_yes <- my.frame[my.frame$target == "YES", ]
T_no <- my.frame[my.frame$target == "NO", ]

# Histogramas
# Slope
slope_yes <- hist(T_yes$slope, plot=FALSE)
slope_no <- hist(T_no$slope, plot=FALSE)

plot(slope_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Slope", ylab="Frecuencia")
plot(slope_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topleft", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# CA
ca_yes <- hist(T_yes$ca, plot=FALSE)
ca_no <- hist(T_no$ca, plot=FALSE)

plot(ca_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="CA", ylab="Frecuencia")
plot(ca_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topright", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# Thal
thal_yes <- hist(T_yes$thal, plot=FALSE)
thal_no <- hist(T_no$thal, plot=FALSE)

plot(thal_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Thalassemia", ylab="Frecuencia")
plot(thal_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topleft", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# Age
age_yes <- hist(T_yes$age, plot=FALSE)
age_no <- hist(T_no$age, plot=FALSE)

plot(age_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Age", ylab="Frecuencia")
plot(age_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topright", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# Trestbps
trestbps_yes <- hist(T_yes$trestbps, plot=FALSE)
trestbps_no <- hist(T_no$trestbps, plot=FALSE)

plot(trestbps_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Trestbps", ylab="Frecuencia")
plot(trestbps_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topright", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# Chol
chol_yes <- hist(T_yes$chol, plot=FALSE)
chol_no <- hist(T_no$chol, plot=FALSE)

plot(chol_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Chol", ylab="Frecuencia")
plot(chol_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topright", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# Thalach
thalach_yes <- hist(T_yes$thalach, plot=FALSE)
thalach_no <- hist(T_no$thalach, plot=FALSE)

plot(thalach_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Thalach", ylab="Frecuencia")
plot(thalach_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topright", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

# Oldpeak
oldpeak_yes <- hist(T_yes$oldpeak, plot=FALSE)
oldpeak_no <- hist(T_no$oldpeak, plot=FALSE)

plot(oldpeak_yes, col=rgb(0, 0, 1, 1/4), main=NULL, xlab="Oldpeak", ylab="Frecuencia")
plot(oldpeak_no, col=rgb(1, 0, 0, 1/4), add=T)
legend(x = "topright", legend=c("Yes", "No"),  
       fill = c(rgb(0, 0, 1, 1/4), rgb(1, 0, 0, 1/4)))

Para el caso de las variables ca, thal y slope que no tienen valores dento de un rango en sí, se tiene que no se observa muy bien las distribuciones a comparación de los otros gráficos, sin embargo, igualmente es apreciable.

Pregunta 2: Teorema Central del Limite

Pruebe el teorema central del limite aplicando un muestreo de la media en las distribuciones Gamma, Normal y una a su elección. Grafique los resultados obtenidos y señale aproximadamente el numero de muestreos necesarios para obtener el resultado esperado, pruebe esto con las siguientes cantidades de muestreo \(\{10,100,1000,5000\}\). ¿El efecto ocurre con el mismo número de muestreo para todas las distribuciones?.

Respuesta

# Definicion de variables o estructuras necesarias para el muestreo.
dist1 <-  rgamma(10000, shape=2)
dist2 <- rnorm(10000)
dist3 <- rbinom(10000, size=1, prob=3/5)

# Realizar el muestreo de las distribuciones.
gen.muestra <- function(m){
  result <- list(c(),c(),c())
  for(i in 1:m) {
    n <- 1000
    muestra1<-sample(dist1,n)
    mean1 <- mean(muestra1)
    muestra2 <- sample(dist2,n)
    mean2 <- mean(muestra2)
    muestra3 <- sample(dist3,n)
    mean3 <- mean(muestra3)
    result <- list(c(result[[1]],mean1),c(result[[2]],mean2),c(result[[3]],mean3))
  }
  return(result)
}
for (i in c(10, 100, 1000, 5000)) {
  resultado.exp <- gen.muestra(i)
  hist(resultado.exp[[1]],
     main = paste("Distribución gamma con", i, "muestras"),
     xlab = "valores",     
     ylab = "frecuencia",  
     breaks = seq(min(resultado.exp[[1]]), max(resultado.exp[[1]]), length.out = 15)
  )
  hist(resultado.exp[[2]],
     main = paste("Distribución normal con", i, "muestras"),
     xlab = "valores",     
     ylab = "frecuencia",
     breaks = seq(min(resultado.exp[[2]]), max(resultado.exp[[2]]), length.out = 15)
  )
  hist(resultado.exp[[3]],
     main = paste("Distribución binomial con", i, "muestras"),
     xlab = "valores",     
     ylab = "frecuencia",
     breaks = seq(min(resultado.exp[[3]]), max(resultado.exp[[3]]), length.out = 15)
  )
}

La tercera distribución seleccionada fue la binomial.

Según los resultados obtenidos en los gráficos, el valor que nos entrega lo esperado del teorema sería 5000. Si bien con mil ya se obtienen distribuciones muy cercanas a la normal, aún se notan unas pequeñas diferencias, las cuales ya casi desaparecen con el valor 5000.

Con respecto a la pregunta, el efecto ocurre con todas las distribuciones, siempre y cuando sea un valor adecuado para ello (por ejemplo, con 10, en los gráficos no se aprecia ninguna distribución normal).

Pregunta 3: Ley de los Grandes Numeros.

Lanzamiento de monedas

Realice el experimento de lanzar una moneda cargada 1000 veces y observe el comportamiento que tiene la probabilidad de salir cara. Para realizar el experimento considere que la moneda tiene una probabilidad de \(5/8\) de salir cara. Grafique el experimento para las secuencias de intentos que van desde 1 a 1000, señalando el valor en que converge la probabilidad de salir cara.

Respuesta

prob <- 5/8
result <- c()

# Simular lanzamientos, 1 = cara
for (lanzamientos in 1:1000) {
  num <- rbinom(1, 1, prob)
  result <- c(result, num)
}

# Primero calculamos la cantidad de 1s en cada punto del vector de resultados
cant.unos <- cumsum(result)
print(result[0:10])
##  [1] 0 0 1 1 1 1 0 0 1 0
print(cant.unos[0:10])
##  [1] 0 0 1 2 3 4 4 4 5 5
# Después calculamos la proporción de 1s en cada punto del vector de resultados
# (normalizamos)
prop.unos <- cant.unos/seq(from=1, to=length(cant.unos), by=1)
print(prop.unos[0:10])
##  [1] 0.0000000 0.0000000 0.3333333 0.5000000 0.6000000 0.6666667 0.5714286
##  [8] 0.5000000 0.5555556 0.5000000
# Gráfico de la convergencia
plot(prop.unos, xlab="Lanzamientos", ylab="Probabilidad")
abline(h=5/8, col="red")

El valor en el que converge la probabilidad de salir cara es a los 200 lanzamientos aproximandamente. Esto, ya que, antes de este valor se ve un comportamiento caótico, el cual luego converge a 0.625 a partir del lanzamiento indicado.

El problema de Monty Hall

Remontándonos en la televisión del año 1963, en USA existía un programa de concursos donde los participantes debían escoger entre 3 puertas para ganar un premio soñado. El problema del concurso era que solo detrás de 1 puerta estaba el premio mayor, mientras que detrás de las otras dos habían cabras como “premio”.

Una de las particularidades de este concurso, es que cuando el participante escogía una puerta, el animador del show abría una de las puertas que no fue escogida por el participante (Obviamente la puerta abierta por el animador no contenía el premio). Tras abrir la puerta, el animador consultaba al participante si su elección era definitiva, o si deseaba cambiar la puerta escogida por la otra puerta cerrada. Un vídeo que puede ayudar a comprender mejor este problema en el siguiente link.

Imagine que usted es participante del concurso y desea calcular la probabilidad de ganar el gran premio si cambia de puerta en el momento que el animador se lo ofrece. Utilizando listas/arrays/vectores simule las puertas del concurso, dejando aleatoriamente el premio en alguna posición del array. Hecho esto, genere un numero de forma aleatoria para escoger una de las puerta (posiciones de la estructura), para luego ver si cambiando de posición tendrá mayores posibilidades de ganar el premio. Genere N veces el experimento y grafique cada una de las iteraciones, tal como se hizo en el ejercicio de las monedas.

Respuesta:

# Creamos una función que simule el juego
montyhall <- function(cambiar = TRUE) {
  puertas <- sample(1:3, 3)  # Puertas donde la posición que tiene el 3 es el premio
  posicion <- sample(1:3, 1) # Elección del participante.
  if (cambiar) {
    eleccion <- puertas[posicion] != 3
  } else {
    eleccion <- puertas[posicion] == 3
  }
  ifelse(eleccion, return(1), return(0)) # Retornamos la elección, esta puede que tenga el premio o no
}

# Función que simula N juegos
n_juegos <- function(n = 1000, cambiar_puerta = TRUE){
  results <- c()
  while (n > 0) {
    result.game <- montyhall(cambiar_puerta)
    results <- c(results, result.game)
    n <- n-1
  }
  return(results)
}

plot(cumsum(n_juegos())/seq(from=1, to=length(n_juegos()), by=1), xlab="Juegos", ylab="Probabilidad", main="Experimento de Monty Hall al cambiar de puerta")
abline(h=2/3, col="red")

plot(cumsum(n_juegos(cambiar_puerta=F))/seq(from=1, to=length(n_juegos()), by=1), xlab="Juegos", ylab="Probabilidad", main="Experimento de Monty Hall al mantener la puerta")
abline(h=1/3, col="red")

Según el gráfico generado, se concluye que al cambiar la puerta efectivamente aumenta la probabilidad de ganar al doble a que si mantuviese su elección inicial. Esto ya que la probabilidad llega a ser aproximadamente 2/3 y en el otro caso, esta es de solo 1/3.

Pregunta 4: ¿Independencia?

Ustedes disponen de los dados D1 y D2, los cuales no están cargados y son utilizados para comprobar que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\) cuando el evento A es independiente del B. Para estudiar la independencia considere que los eventos A y B se definen de la siguiente manera; sea A el evento dado por los valores obtenidos en el lanzamiento del dado D1, este está compuesto por \(A=\{D1=1,D1=2,D1=6\}\). Por otro lado, el evento B viene dado por los valores obtenidos con el dado D2, el que está conformado por \(B=\{D2=1,D2=2,D2=3,D2=4\}\). Con esto, tendremos un \(\mathbb{P}(A)=1/2\), \(\mathbb{P}(𝐵)=2/3\) y \(\mathbb{P}(AB)=1/3\). Compruebe de forma gráfica que al realizar 1000 lanzamientos (u otro valor grande que usted desea probar) se visualiza que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\).

Hecho lo anterior, compruebe el comportamiento de un segundo grupo de eventos, dados por el lanzamiento de solo el dado D1. Donde, los eventos para D1 quedan definidos como: \(A =\{D1=1,D1=2,D1=6\}\) y \(B=\{D1=1,D1=2,D1=3\}\). ¿Se observa independencia en este experimento?.

Se le aconseja que para simular los lanzamientos de dados utilice la función sample() para generar valores aleatorios entre 1 y 6. Compruebe los números generados por la función con los casos favorables de cada uno de los eventos a ser estudiados.

Respuesta:

# Primer grupo de eventos
N_lan = 1000

L_A = c(1, 2, 6)
L_B = c(1, 2, 3, 4)

v_A = sample(1:6, N_lan, TRUE)
v_B = sample(1:6, N_lan, TRUE)

event <- function(v, l) {
  ifelse(v %in% l, return(1), return(0))
}

occ_A <- sapply(v_A, event, L_A)
occ_B <- sapply(v_B, event, L_B)
occ_AB <- occ_A & occ_B

prob_A <- cumsum(occ_A) / seq(1:N_lan)
prob_B <- cumsum(occ_B) / seq(1:N_lan)
prob_AB <- cumsum(occ_AB) / seq(1:N_lan)

plot(prob_A, col="red", xlab="Lanzamientos", ylab="Probabilidad", main="Independencia de eventos", ylim=c(0, 1))
lines(prob_B, col="blue", type="p")
lines(prob_AB, col="green", type="p")
legend("topright", c("Probabilidad A", "Probabilidad B", "Probabilidad AB"), lty=1:3, col=c("red", "blue", "green"), cex=.80)

Para el primer experimento es posible evidenciar como el evento A y B convergen a sus probabilidades teóricas, que en este caso son 1/2 y 2/3 respectivamente, y además se observa como la probabilidad de que ocurran ambos eventos simultáneamente también tiende al valor esperado dado que efectivamente son eventos independientes, que en este caso corresponde a \(\mathbb{P}(AB) = \mathbb{P}(A) \cdot \mathbb{P}(B) = \frac{1}{2} \cdot \frac{2}{3} = \frac{1}{3}\).

# Segundo grupo de eventos
N_lan = 1000
  
L_A = c(1, 2, 6)
L_B = c(1, 2, 3)

v_A = sample(1:6, N_lan, TRUE)

occ_A <- sapply(v_A, event, L_A)
occ_B <- sapply(v_A, event, L_B)
occ_AB <- occ_A & occ_B

prob_A <- cumsum(occ_A) / seq(1:N_lan)
prob_B <- cumsum(occ_B) / seq(1:N_lan)
prob_AB <- cumsum(occ_AB) / seq(1:N_lan)

plot(prob_A, col="red", xlab="Lanzamientos", ylab="Probabilidad", main="Independencia de eventos", ylim=c(0, 1))
lines(prob_B, col="blue", type="p")
lines(prob_AB, col="green", type="p")
legend("topright", c("Probabilidad A", "Probabilidad B", "Probabilidad AB"), lty=1:3, col=c("red", "blue", "green"), cex=.80)

Para el segundo experimento en cambio, es posible evidenciar que no son independientes puesto que la probabilidad a la que converge la intersección de estos eventos es distinta a la esperada teóricamente en caso de que lo fueran, que en este caso corresponde a 1/4, lo cual es distinto al valor que se obtiene al calcular esta probabilidad considerando que son dependientes, la cual es \(\mathbb{P}(AB) = \mathbb{P}(A) + \mathbb{P}(B) - \mathbb{P}(A \vee B) = \frac{1}{2} + \frac{1}{2} - \frac{2}{3} = \frac{1}{3}\), que efectivamente es el valor al que se acerca este experimento.

Pregunta 5: La Ruina del Jugador

Un amigo ludópata suyo le comenta que el truco de jugar en el casino está en no parar de apostar y apostando lo mínimo posible. Ya que así, tienes mas probabilidades de ganar el gran pozo que acumula el juego. Usted sabiendo la condición de su amigo, decide no creer en su conjetura y decide probar esto a través de un experimento.

Para realizar el experimento usted decide asumir las siguientes declaraciones, bajo sus observaciones:

En el primer experimento deberá obtener la evolución de los fondos hasta que el jugador se queda sin fondos para jugar. Puede ser útil seguir la lógica de una moneda cargada para realizar esto. Pruebe esto con una apuesta igual a 5, 25 y 50 graficando los resultados. Comente los resultados obtenidos.

Para la segunda parte de este experimento, con las funciones ya creadas, proyecte 5000 apuestas y obtenga la probabilidad a la que converge el experimento empezando con una apuesta de: 5, 25 y 50. Para este experimento considerar como éxito (1) cuando su función ruina supera los 200 y considere lo contrario cuando su función perdida es menor o igual a 0.

Respuesta

# Función para obtener el desarrollo de las apuestas
ruina <- function(fondos = 100, apuesta = 5){
  win_prob <- 8/19
  vec_fondos <- c()
  while (0<fondos & fondos<200) {
    value <- runif(1, 0, 1)
    if (value <= win_prob) {
      fondos <- fondos + apuesta
    } else {
      fondos <- fondos - apuesta
    }
    vec_fondos <- append(vec_fondos, fondos)
  }
  return(vec_fondos)
}

plot(ruina(), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 5)")

plot(ruina(apuesta = 25), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 25)")

plot(ruina(apuesta = 50), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 50)")

Con estos resultados empíricos es posible evidenciar que la estrategia propuesta por el amigo no es factible puesto que finalmente termina perdiendo sus fondos, esto debido a que la probabilidad de ganar sea menor a un medio, generando que la tendencia vaya a perder dinero. La principal diferencia en la cantidad de dinero apostada influye en que decrece la cantidad de juegos a realizar, ya que aumenta o disminuye más rápidamente el dinero entre apuestas.

tries <- 5000

results_5 <- c()
results_25 <- c()
results_50 <- c()
for (i in 1:tries) {
  results_5 <- append(results_5, ifelse(tail(ruina(apuesta=5), n=1) <= 0, 0, 1))
  results_25 <- append(results_25, ifelse(tail(ruina(apuesta=25), n=1) <= 0, 0, 1))
  results_50 <- append(results_50, ifelse(tail(ruina(apuesta=50), n=1) <= 0, 0, 1))
}

prob_5 <- cumsum(results_5) / seq(1:tries)
prob_25 <- cumsum(results_25) / seq(1:tries)
prob_50 <- cumsum(results_50) / seq(1:tries)

plot(prob_5, xlab="Cantidad de juegos", ylab="Probabilidad")

plot(prob_25, xlab="Cantidad de juegos", ylab="Probabilidad")

plot(prob_50, xlab="Cantidad de juegos", ylab="Probabilidad")

Finalmente, con estos gráficos es posible evidenciar que el consejo del amigo estaba errado en 2 cosas, la primera es que la tendencia está claramente en perder contra el casino, y la segunda es que conviene más apostar una mayor cantidad de dinero.

A work by CC6104